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ABSTRACT 

We have compared molecular line emission to dust continuum emission and modeled molecular lines 
using Monte Carlo simulations in order to study the depletion of molecules and the ionization fraction 
in three preprotostellar cores, L1512, L1544, and L1689B. L1512 is much less dense than L1544 and 
L1689B, which have similar density structures. L1689B has a different environment from those of L1512 
and L1544. We used density and temperature profiles, calculated by modeling dust continuum emission 
in the submillimeter, for modeling molecular line profiles. In addition, we have used molecular line 
profiles and maps observed in several different molecules toward the three cores. We find a considerable 
diversity in chemical state among the three cores. The molecules include those sensitive to different 
timescales of chemical evolution such as CCS, the isotopes of CO and HCO+, DCO"*", and N2H"'". The 
CO molecule is significantly depleted in L1512 and L1544, but not in L1689B. CCS may be in the 
second enhancement of its abundance in L1512 and L1544 because of the significant depletion of CO 
molecules. N2H^ might already start to be depleted in L1512, but it traces very well the distribution 
of dust emission in L1544. On the other hand, L1689B may be so young that N2H+ has not reached its 
maximum yet. The ionization fraction has been calculated using H^^CO"^ and DCO^. The result shows 
that the ionization fraction is similar toward the centers of the three cores. This study suggests that 
chemical evolution depends on the absolute timescale during which a core stays in a given environment 
as well as its density structure. 



1. INTRODUCTION 

Many studies have been done on low mass-star forma- 
tion in detail, theoretically and observationally, over the 
past two decades. These studies were possible because 
many low mass protostellar cores are close, isolated, and 
geometrically simple. In addition, they evolve slowly, al- 
lowing the process to be observed in distinct stages. The 
information accumulated by those studies now provides 
the big picture of the evolutionary stages of low-mass star 
formation. Lada (1987), Myers & Ladd (1993) and Andre 
et al. (1993) developed an evolutionary sequence based 
on the spectral energy distribution (SED) of dust contin- 
uum emission through Class to Class III sources, and 
Adams et al. (1987) supported the evolutionary sequence 
by modeling the SED theoretically. In addition, several 
groups (Larson 1969, Penston 1969, Shu 1977, Shu et 
al. 1987, Foster & Chevalier 1993, Ciolek & Mouschovias 
1994, McLaughlin & Pudritz 1997) have developed dy- 
namical evolution models that theoretically predict differ- 
ent structures of density and velocity in protostellar cores, 
based on differing initial conditions. 

In spite of the big picture of low-mass star formation, 
the initial collapse is poorly understood. The collapse dy- 
namics and the timescale of evolution are crucially depen- 
dent on initial conditions. Preprotostellar cores (hereafter 
PPCs) are believed to be gravitationally bound, but they 
have no central hydrostatic protostar (Ward-Thompson et 
al. 1994, 1999). Therefore, PPCs can be considered as 



the potential sites of future star formation and give us 
chances to probe the initial conditions of star formation. 
The series of papers of Ward-Thompson et al. (1994, 1999, 
2000, and 2001) and the paper of Andre et al. (1996) ar- 
gued that the density of PPCs is relatively constant in the 
inner cores and falls off at larger radii. This is very differ- 
ent from the singular isothermal sphere that leads to the 
"inside-out" collapse (Shu 1977). However, they assumed 
a uniform temperature in PPCs to calculate density struc- 
tures. Evans et al. (2001) and Zucconi et al. (2001) show 
that the temperature of a dust grain decreases toward the 
centers of PPCs, assuming that dust grains are heated 
only by the attenuated external radiation field and that 
the radiation from the core itself is optically thin. Based 
on their results, Evans et al. (2001) argue that the con- 
tinuum emission is significantly dependent on the temper- 
ature structure as well as the density structure, and they 
show that Bonnor-Ebert spheres can fit well the submil- 
limeter dust continuum emission in PPCs. In agreement, 
the recent paper of Ward-Thompson et al. (2002) has 
shown that the color temperatures of PPCs that are cal- 
culated from the 170 and 200 fj,m dust emission observed 
by ISOPHOT decrease toward their centers. As a result 
of the lower temperature at the center, the density at the 
center must be more peaked than is the intensity of dust 
continuum emission. Evans et al. (2001) showed that even 
a power law could fit the density structure of L1544. 
Ciolek & Mouschovias (1994) modeled ambipolar diffu- 
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sion for the quasi-static evolution of a uniformly magne- 
tized isothermal molecular cloud, finding relatively con- 
stant density toward the center. However, there has 
been debate about the importance of ambipolar diffusion 
because of differences between inward velocity observed 
(Tafalla et al. 1998, WiUiams et al. 1999) and calcu- 
lated from the magnetic field (Goodman and Heiles 1994, 
Andre et al. 1996, Ward-Thompson et al. 1999, Crutcher 
and Troland 2000, Ward- Thompson et al. 2000). Those 
differences might be resolved by calciilating more accurate 
density structure and, in turn, the correct temperature 
structure from the observed data. Ciolek and Basu (2000) 
showed if the density at the center of a PPC were higher, 
then their model could explain the observed magnetic field 
and inward velocity. 

In addition to density and temperature distributions, we 
need to understand the ionization structure because mag- 
netic relaxation by ambipolar diffusion is tightly related 
to the ionization structure (Rawlings 2000). Ionization 
fraction can affect the process of gravitational collapse be- 
cause the coupling between ions and magnetic field plays a 
key role in preventing free-fall collapse of a dense core and 
slowing the drift of neutral gas through the magnetic field 
(ambipolar diffusion). In addition, the gas-phase chem- 
ical reactions depend significantly on ionization fraction, 
and the ionization fraction is highly sensitive to the de- 
pletion of molecules because some molecules play a cru- 
cial role in destroying ions. The degree of ionization in 
cold dense cores can be measured by the [DCO+]/[HCO+] 
abundance ratio (Watson 1977), which is directly related 
to the electron density and depletion fraction. Recently, 
many molecular line observations show evidence for chem- 
ical differentiation and the depletion of some molecules 
such as C^^O, C^^O, and CCS in dense cores (Lada et 
al. 1994, Kuiper et al. 1996, Black et al. 1995, Mc- 
MuUin et al. 1994, Benson and Myers 1998, Willacy et al. 
1998, Kramer et al. 1999, Ohashi et al. 1999, CaseUi et 
al. 1999, Jessop and Ward-Thompson 2001, Tafalla et al. 
2002). They suggested that this trend is due to the sub- 
stantial depletion of molecules onto dust grain surfaces. In 
addition, several groups (Bergin and Langcr 1997, Aikawa 
et al. 2001, Caselh et al. 2002, Li et al. 2002) modeled 
chemical evolution, including gas-dust interaction as well 
as gas-gas interaction, to show the changing depletions of 
some molecules and ionization fraction with time or den- 
sity, and most of them compared their results to real ob- 
servations. Therefore, both gas phase chemical processes 
and dust-gas interactions should be considered in studying 
star formation. 

Those theoretical and observational results suggest that 
the relative distributions of several molecules can be used 
to probe the timescale of the dynamical evolution of a 
PPC. However, Rawlings et al. (1992) and Rawlings and 
Yates (2001) modeled several molecular lines combining 
self-consistently dynamical, chemical, and radiative trans- 
fer models, and showed that the abundance variation and 
the line profiles are very sensitive to the free parameters in 
the chemical model. Consequently, we need to use multiple 
transitions and maps in order to constrain these models. 

In this paper, we explore the chemical and dynamical 
differences in three PRCs (L1512, L1544, and L1689B) 
comparing molecular depletion and ionization fraction by 



simple analysis of integrated intensity and column density 
and by modeling several molecular line profiles in detail. 
In this study, we cover a very early-time molecule (CCS), 
an early-time molecule (CO), a middle-time molecule 
(HCO^), and a very late-time molecule (N2H+) in chemi- 
cal evolution. According to the chemical models of Bergin 
and Langer (1997) and Aikawa et al. (2001) that used an 
initial density of 10'' cm~^, the CCS abundance decreases 
rapidly after a few times 10^ years since CCS is a polar 
molecule, and it is easily frozen out onto the surface of dust 
grains. In addition, CCS is destroyed by HCO"*" or H3O+, 
whose abundances increase until about 10^ years. Even 
though CO is a non-polar molecule, it has very high bind- 
ing energy (1740 K onto water ice mantle grains, Bergin 
and Langer 1997). As a result, CO starts to be depleted 
onto the surface of dust grains after about 10^ years caus- 
ing the depletion of HCO+ in turn. HCO+ is formed by 
the interaction between and CO or H2D+ and CO. 
However, N2, which is the precursor of N2H+, starts to be 
abundant much later than other molecules and reaches its 
maximum around 10^ years. N2 has low binding energy 
(750 K, Aikawa et al. 2001) onto the surface of dust grains, 
so it is not frozen out from the gas phase easily. In addi- 
tion, H3 , which is the other precursor of N2H+, increases 
as CO is depleted. Another cause of the enhancement of 
N2H"'" at late times is the decline of the destruction of 
N2H+ by CO, which is significantly depleted at late times. 

We observcid multiple positions in each core to study the 
chemical distribution within a core in addition to study- 
ing the chemical differences among the three PRCs. The 
most important aspect of this study is that we start with 
reasonable physical models calculated from dust analysis. 
The density and temperature structures found by Evans et 
al. (2001) are adopted for modeling molecular line profiles. 

In section 2, we summarize observational parameters, 
and section 3 shows the resiilts of observations. In sec- 
tion 4, we analyze column densities of several molecules. 
Section 5 compares modeled molecular line profiles and 
observed line profiles. Finally, we discuss our results re- 
garding depletion and ionization in the view of chemical 
and dynamical evolution in section 6. 

2. OBSERVATION 

We observed the three pre-protostellar cores in C^^O 

(J = 2 - 1 and J = 3 - 2), C^^O (J = 2 - 1), DCO+(J = 
3 - 2), HCO+ (J = 3 - 2), and Hi^^CO+ (J = 3 - 2) with 
the 10.4 m telescope of the Caltech Submillimeter Obser- 
vatory (CSO) at Mauna Kea, Hawau from 1995 to 2002. 
We used an SIS receiver with an acousto-optic spectrom- 
eter (AOS) with a 50 MHz bandwidth and 1024 channels, 
and the frequency resolution was about 2 to 2.5 channels, 
0.13 to 0.16 km s"^ at around 220 GHz. The pointing un- 
certainty was about 4" on average. We also observed these 
cores in Hi3C0+ (J = 1 - 0), N2H+ (J = 1 - 0), and CCS 
(Nj =43 — 82) with the 45 m telescope of the Nobeyama 
Radio Observatory in Japan on January 2002. We used the 
S40, S80, and SlOO receivers with an AOS with 40 MHz 
bandwidth and 2048 channels. The frequency resolution 
was about 2 channels (about 40 kHz) so that the velocity 
resolution was about 0.27, 0.14, and 0.13 km at 43, 
86, and 100 GHz, respectively. The observed sources are 
listed in Table 1; the observed lines and the main beam 
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efBciency {rjmb) of each observation are summarized in Ta- 
ble 2. In Table 1, the last column shows the offsets of the 
dust peaks from the given coordinates; all offsets in figures 
other than Figure 2 and 3 are relative to the dust peak, 
assumed to be the center of the core in all modeling. 

3. RESULTS 

3.1. Spectra 

Figure 1 shows the spectra of molcciilar lines observed 
toward the centers of the three PPCs. The centers are the 
positions of the centroid of the dust emission (Shirley et al. 
2000). Most of the molecular lines have not been mapped 
completely but observed through two orthogonal cuts in 
a — (5, crossing at the centers of dust emission in the cores. 
We do not see any consistent trend in the intensity of each 
spectrum from core to core. For example, C^^O and C^^O 
lines are the strongest in L1689B, but N2H+ and CCS are 
the strongest in L1544. All the lines except for N2H+ and 
CCS are much weaker in L1512 than in L1689B. The line 
width (FWHM) of the C^^O J = 2-1 line is about 0.3, 0.5, 
and 0.6 km s'^ in L1512, L1544, and L1689B, respectively. 

3.2. Maps 

Figure 2 shows maps of three PPCs, where the grey 
scale indicates 850 /im dust continuum intensity and the 
contours are for integrated intensities of C^^O J = 2 — 1 
and DCO"'" J = 3 — 2 lines. Wc can see that the peaks 
of dust emission and C^^0J = 2 — 1 emission are not 
consistent. In the map for L1512, especially, the centroid 
of the dust emission coincides with a hole in the C'^^O 
J = 2 — 1 emission. In contrast, the distribution of DCO"*" 
J = 3 — 2 emission appears to peak toward the center of 
dust emission in L1544 and L1689B, even though our maps 
of 000+ J = 3 — 2 emission are not complete. 

4. SIMPLE ANALYSIS 

In this section, we use a simple analysis with standard 
techniques for comparison to similar work done by others, 
and, in the next section, we will deal with detailed models. 
We use dust emission as the tracer of column density least 
affected by chemistry in this section, and we assume that 
the dust emissivities do not change with radius. If tem- 
perature decreases toward the center, the column density 
calculated from dust emission with the assumption of an 
isothermal case could be underestimated toward the cen- 
ter. Even though the absolute column density calculated 
from the dust emission can change with the assumed dust 
opacity by a factor of 2 or 3, its relative distribution with 
radius can be traced relatively well by dust emission. 

4.1. Integrated Intensities 

We can compare the distribution of a molecular com- 
ponent with that of the dust component within a core 

using the integrated intensity of an optically thin molec- 
ular line and the dust continuum intensity observed at a 
submillimeter wavelength. In the upper panels of Figure 
3, wc plot the logarithmic integrated intensities of CCS 
Nj = 43 - 32, H13C0+ J = 1 - 0, and N2H+ J = 1 - lines 
in order to compare to the 850 nm dust continuum inten- 
sity through the cuts with constant declination shown in 
Figure 2. The dust continuum intensity and the H^^CO+ 



intensity have been shifted upward to compare the shapes 
of intensity distributions with radius more effectively, and 
the intensity of dust continuum has been measured into 
apertures of 20", spaced by 20" for consistency with the 
Hi^''CO+ and N2II+ observations. 

The dust emission of L1512 is about half of that of 
LI 544, and the distribution is not as peaked as that of 
L1544 and L1689B. L1544 and L1689B have similar dust 
emission even though it is more peaked in LI 544. How- 
ever, the CCS intensity of L1689B is more similar to that 
of L1512 rather than to that of L1544, and the CCS inten- 
sity of L1544 is about twice those in L1512 and L1689B. 
The CCS intensity in all the three PPCs is not peaked as 
much toward the center as the dust emission. However, it 
has a bump toward the centers of L1512 and L1544 (see 
§6.1). In L1689B, the bump is not certain because of the 
low signal-to-noise ratio. 

If we check the relative variations of the integrated in- 
tensity of the H^^CO+ J = 1 — line with radius in each 
core, in L1689B, it is skewed to the west but follows well 
the dust continuum emission, decreasing at the edge of 
the core. On the other hand, the distribution of the inte- 
grated intensity of the Hi'^CO+ J = 1 - hne in L1512 is 
not well matched with dust continuum, and the integrated 
intensity is even depressed toward the center. We also see 
the depression of H^^CO"*" emission at the center of L1544. 
This depression was predicted by Caselli et al. (2002) in 
their chemical models, but not observed, probably because 
their angular resolution (28") was worse than ours (18"). 
We suspect that the optical depth of the H^^CO+ J = 1-0 
line in L1544 and L1689B is not negligible. Therefore, the 
depression in those cores could be partially caused by op- 
tical depth. However, L1689B does not show a depression 
even though it has a density structure similar to LI 544 
suggesting that depletion is playing a major role in L1544. 

The N2H"'" intensity is well peaked in L1544 and shows 
almost the same shape as the dust emission. L1689B also 
shows similar shapes in dust emission and N2H+ intensity. 
However, the N2H+ intensity of L1689B is only two thirds 
of that of LI 544, even though these two cores have very 
similar dust emission and density structures, with central 
densities of 10® cm""^. In L1689B, the peak of the intensity 
of N2H+ is shifted from that of the intensity of the dust 
emission similarly to the intensity of H^'^CO+. However, 
the ratio of the intensity of H^^CO"'" to the intensity of 
N2H+ at the center of L1689B, is greater than that in the 
other two cores. 

In L1512, with a central density of 10^ cm'^, N2H+ 
emission is more extended than the dust, unlike in the 
other cores. Since the optical depth of the N2H+ J = 1 — 
line in L1512 is smaller than in L1544 (Benson et al. 1998), 
this might indicate the possible presence of depletion at 
the dust peak. 

In the following section, we compare the distributions of 
CO isotopes to that of dust emission by using molecular 
hydrogen column densities calculated from them instead of 
using the intensity distribution because the optical depth 
effect in the C^^O line can be calculated, and the standard 
abundance of CO is known. 



4.2. Column Density 
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We compare the H2 column densities calculated from ob- 
served molecular line intensities with that calculated from 
submillimeter continuum emission. This simple analysis 
assumes that 1) dust continuum emission traces all mate- 
rial along a line of sight in a core, and dust and gas arc 
well mixed, 2) the core is isothermal, and dust and gas 
are well coupled so that they have same temperature, 3) 
molecular abundances are constant through the core, and 
4) all levels are in LTE. Therefore, the excitation temper- 
atures in all levels arc the same as the kinetic temperature 
and the same as the dust temperature. 

4.2.1. The H2 column density from molecular emission 

lines 

If we assume a line is very optically thin so that absorp- 
tion can be completely neglected, the equation of radiative 
transfer is given by 



ds 



An 



-rin 



(1) 



where n„ is the number density of upper level, and A^i is 
the Einstein coefficient. For the transition J ^ J — 1 in a 
linear molecule, the resulting relation between N{x), the 
column density of molecule x, and the integrated intensity 
of the line, is 



N{x) = 



Tr dv, 



(2) 



where £■„ = hBJ{J + 1), B is the rotational constant, 
Tex is the excitation temperature above the ground state, 
fj, is the dipole moment, and Q is the partition function. 
Equation 2 is valid only in the limit that t ^ 0. For finite 
optical depth, an optical depth correction can be applied: 
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In principle, one should also correct for the presence of 
the cosmic background radiation, Tcmb, but for the lines 
considered here, at A < 1.4 mm, the Planck function in 
temperature units is less than 0.23 K, completely negligi- 
ble compared to likely values oiT^x- The effect on N{x) of 
neglecting Tcmb is less than 5% for Tg^ > 10 K and r < 1, 
typical of the conditions relevant here. These effects are 
less than any plausible calibration uncertainty. 

The rotational constant and the dipole moment of each 
molecule used for this calculation are listed in Table 3. 
For the excitation temperature, we used the kinetic tem- 
perature of 10 K, a typical value in dark molecular clouds, 
under the assumption that all lines are thermalizcd, that 
is, T, = Tex = Tfe = 10 K. If we use 20 K for this calcula- 
tion, the result is changed only by a factor of 1.2 for the 
J = 2 - 1 line of Ci^O and C^^O. 

We calculated the column densities using C^^O J = 2 — 1 
and C-'^^O J = 2 — 1 lines. First, the lines were assumed to 
be optically thin. However, we can calculate the optical 
depth of the J = 2 - 1 line directly by comparing the 
relative intensities of its 9 hyperfine components (Ladd et 
al. 1998). The optical depth of the C"0 J 2 - 1 Hue 
was obtained by a fitting method in the program CLASS. 
The panels in the third row of Figure 1 shows the results 
of fitting the hyperfine components of the J = 2 — 1 



line in the center of each core. The total optical depths 
for L1544 and L1689B are not small enough to be ignored 
(Table 4), even though each individual hyperfine compo- 
nent is optically thin. In order to correct for this optical 
depth effect, wc applied equation (3) to the above calcu- 
lation for an optically thin line. For the C^^O J = 2 — 1 
line, Tv is the optical depth of the line center that is cal- 
culated by multiplying the relative strength (0.693) of the 
main group, which has four close components, by the to- 
tal optical depth. The correction factors at the centers of 
L1544 and L1689B are about 1.4. Note that we do not 
use the excitation temperature of C^''0 J = 2 — 1 obtained 
by fitting the hyperfine structure to calculate the column 
density of C^''0 because this method H-ssiiinGS t licit Tqj; is 
constant along the line of sight. Our detailed models (see 
§5.3.1) shows that this assumption leads to incorrect re- 
sults. Also we want to use the same method for C^^O and 
C-'^^O and compare the results. 

Finally, we convert the column density of each molecule 
into the H2 column density using the abundance {X) of 
molecule x; 



iV(H2) 



N{x) 



(4) 



We used a CO abundance of 2.7 x 10"'' (Lacy et al. 1994), 
about three times bigger than the abundance used in other 
studies such as Bacmann et al. (2002) and Jorgensen 
et al. (2002). The abundance ratios [CO/C^^O] and 
piSQ/ci^o] are about 560 (± 25) and 3.2 (± 0.2), re- 
spectively (Wilson and Rood, 1994). The abundance of 
each molecule used in calculating the H2 column density 
is listed in Table 3. 

4.2.2. The H2 column density from dust continuum 

emission 

If the emission from dust is optically thin, the observed 
flux density {Si,) can be related to the column density of 

gas by 



N{-R2) = 



(5) 



where fi is the mean molecular weight, mn is the atomic 
mass unit, is the mass opacity of dust per gram of gas, 
Bi, is the Planck function, and f7 is the aperture solid an- 
gle. — {-kO'^) / {Aln2) for a circular Gaussian aperture. 

We use column 5 of the Ossenkopf-Henning dust opacity 
table (Ossenkopf & Henning 1994), which represents ag- 
glomerated dust grains with thin ice mantles, and assume a 
gas-to-dust mass ratio of 100 to get K850 ~ 0.018 cm^ g~^. 
We use the 850 jim data from Shirley et al. (2000) but 
measure the flux in an aperture of 33", the size of main 
beam in the observations of C^^O. If the dust tempera- 
ture. Til is 10 K, the H2 column density can be calculated 
from 



N{H2) = 3.36 X lO-"^ X 5850 (Jy) cm 



(6) 



If the temperature were 20 K, the calculated N(H2) would 
be lower by a factor of 3.3. 
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4.2.3. Results 

C^^O J = 2 - 1 and C^^O J = 2 - 1 lines have been 

used to calculate the H2 column density assuming constant 
molecular abundances. Basically, we assume the lines are 
optically thin. If a molecular line is optically thick, then 
the calculation with equation (2) and (4) underestimates 
the H2 cohimn density. As mentioned in the previous sec- 
tion, however, we can calculate the optical depth of the 
C^^0J = 2— 1 line by fitting its hyperfine structure. We 
have corrected the column density of the C^^O J = 2 — 1 
line for optical depth in L1544 and L1689B even though 
the correction is not very big. 

The results of fitting the hyperfine components of C^^O 
J = 2 — 1 show that the total optical depth of this line 
is about 0.1 in L1512, but it is not negligible at the cen- 
tral positions in L1544 (r w 1.0) and L1689B (t » 0.8). 
In L1512, the error of the calculated optical depth is very 
big, because of the low S /N in the weak components. How- 
ever, the ratio (« 4.3) between the strongest component 
and the second strongest component in L1512 is bigger 
than that 3.3) in the two other cores, so the total opti- 
cal depth would be less than that of L1689B, at least. The 
abundance of C^^O is about 3.2 times greater than that 
of C-'^^O. Therefore, the optical depth of C-'^^O is about 
3.2 greater than the total optical depth of C^^'O, and the 
optical depth of C^'^O .1 = 2 — 1 must be considered in 
calculating column densities. We plot the H2 column den- 
sity calculated from the C^^O J = 2 — 1 line in the lower 
panes of Figure 3, where the optical depth effect has been 
corrected by equation (3). However, in the figure, we use 
the H2 column density calculated from the C"'"*0 J = 2 — 1 
line based on the optically thin line approximation without 
correcting for its optical depth in order to show the signif- 
icance of the optical depth effect. In fact, C-'^^O has been 
observed in fewer positions so that we cannot calculate the 
optical depth of every C^^O spectrum. 

In the figure, we see significant differences between the 
H2 column densities calculated from the 850 /xm emission 
and the C^^O .1 = 2-1 fine in all three PPCs. At cen- 
tral positions, the H2 column density calculated from the 
C^^O J = 2 — 1 line is about 15 times smaller in L1512 and 
L1544, and 8 times smaller in L1689B than the H2 column 
densities calculated from the Ssbo /^ra emission. However, 
after correcting for the optical depth of the C'^^O .1 = 2 — 1 
line, the difference between the H2 column densities cal- 
culated from Ssso and from C^^O has been reduced. In 
L1689B, especially, the difference is just a factor of 2 after 
the correction. 

Another interesting thing is that the distribution of 
C^^O in L1544 does not show the central depression ob- 
served in Caselli et al. (1999). This is probably due to 
the larger beam (33") in this observation compared to 
that (21") of Caselli et al. (1999), similarly to the case 
of H^^CO"*" (see §4.1). This result suggests that inade- 
quate resolution can prevent us for calculating the actual 
depletion of molecules. 

4.3. Conclusions 

The integrated intensities and the H2 column densi- 
ties calculated from several molecules have been compared 
with the intensity of 850 fim dust emission and the H2 
column densities calculated from the dust emission in our 



simple analysis of molecular line observations. The first 
result of this analysis is that the dust emission is more 
peaked than the molecular emission. The depression of 
molecular emission toward the centers is partially caused 
by the optical depth effect on molecular lines. We could 
see this effect comparing the column densities calculated 
from C^^'O J = 2 — 1 (marginally optically thin and cor- 
rected by r) and C^^O J = 2 — 1 (optically thick but not 
corrected by r) lines. The correction for the optical depth 
effect by the C^^O J = 2 — 1 line removes some of the dif- 
ference between the H2 column densities calculated from 
the C^^O J = 2 — 1 line and dust emission. The remaining 
considerable discrepancy would be caused by the depletion 
of CO molecules. This result shows that we can overesti- 
mate the depletion of a molecule if we use an optically thick 
line without correcting for the optical depth. Jorgensen et 
al. (2002) have compared the H2 column densities calcu- 
lated from the 1.3 mm dust contimmm and the C"'^''0 1-0 
line and found similar depletion factors of CO to ours in 
L1544 and L1689B. We will explore the depletion in detail 
through modeling molecular line profiles using the Monte 
Carlo method in the next section. We emphasize that the 
amount of the discrepancy varies from core to core. After 
the correction for the optical depth effect, the H2 column 
density calculated from the C^^O J = 2 — 1 line is 11 in 
L1512, 9 in L1544, and 2 in L1689B times smaller than the 
H2 column density calculated from dust emission (Table 
5). 

The simple method used in this analysis clearly reveals 
evidence for depletion, but it has several limitations for 
quantitative analysis. First, the beam size of each line ob- 
servation is not the same. If we had mapped these cores 
completely in each line, then we would have been able 
to convolve the maps using the same beam size. Second, 
the variations of the excitation temperature and the abun- 
dance of a molcciile along the line of sight arc not included 
in this simple analysis. In addition, the absolute H2 col- 
umn density calculated from dust emission may be uncer- 
tain by a factor of 2 to 3 because of the uncertainty of 
dust opacity. Due to these limitations, we have to model 
molecular lines and compare them to the observed lines in 
order to study depletion more in detail. 

5. DETAILED MODELS 

5.1. Method 

We used the Monte Carlo (hereafter MC) method to 
calculate the radiative transfer of molecular lines (Bernes 
1979, Choi et al. 1995). The MC code for this work has 
been developed by Choi et al. (1995). The MC code 
generates model photons at a random position, in a ran- 
dom direction, and at a random frequency with proper 
random number distributions. These photons go through 
a one-dimensional spherically symmetric molecular cloud 
adjusting the level populations of molecules. The MC code 
calculates the excitation by the model photons and uses 
statistical equilibrium to adjust each level population until 
the criteria of convergence are satisfied. The MC method 
is more powerful than the LVG method because it can deal 
with arbitrary distributions of systematic velocity, density, 
kinetic temperature, microturbulence, and abundance self- 
consistently. Therefore, we can combine reasonable phys- 
ical models from analysis of dust emission with chemical 
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models to produce molecular line simulations, and the de- 
pletion of molecules with radius can be tested. 

For the MC simulation, we needed collision rates at 
lower gas temperatures than are usually available, so we 
extrapolated linearly to 5 K the collision rates for CO 
(Flower and Launay 1985), which we used for C^^O and 
C^^O. Similarly, we extrapolated the downward collision 
rates for HCO+ (Flower 1999) to use for Hi3C0+ and 
DCO"'" . Upward rates were calculated for detailed balance 
for each isotope. The size of each core is fixed at 0.15 pc, 
and microturbulent velocity is assumed constant through 
each modeled core. 

Once the MC code calculates each level population, we 
simulate specific molecular line profiles using the virtual 
telescope (VT) simulation to compare the modeled profiles 
with observed molecular line profiles. We can simulate a 
molecule with multiple transitions and multiple positions 
simultaneously in a specific core with the parameters such 
as beam sizes and main beam efficiencies that we used in 
the actual observations. The remaining lack of realism is 
caused by geometric differences between a 1-dimcnsional 
simulation and the actual core. As seen in Figure 2, the 
three cores are not perfectly spherically symmetric, but 
instead they are elongated. In order to compare simulated 
line profiles with real observations, we averaged the lines 
that have been observed at the same distance from the 
centers, e.g., the lines in the off-positions of ±30" in a 
and S. However, some lines have not been averaged in all 
four directions because the observation did not cover all 
parts. We need to map cores fully and use a 2-dimensional 
code to do more accurate comparisons. 

We model the transitions of four molecules, C^^O (J = 
2-landJ = 3-2),DCO+ (J = 3-2), Hi'''CO+ (J = 1-0 
and J = 3 - 2), and HCO+ (J = 3 - 2). If the gas and 
dust of a core are well mixed and well coupled, we can 
apply the distributions of density and temperature of the 
dust to the gas. However, for a density lower than about 
10^ cm~^, the gas temperature is not coupled well with 
the dust temperature. The excitation temperatures of the 
lines that have high critical density such as 000+ and 
H^"'CO+ are mainly dependent on density so that the in- 
ner, denser region mainly contributes to the intensities of 
those molecular lines, and the outer, less dense region does 
not change our results much. However, we should in prin- 
ciple correct for the difference between gas and dust tem- 
peratures in the outer, less dense parts of cores, especially 
for the molecular lines such as C^^O lines that have low 
critical density. We leave this correction for a future work. 

5.2. Physical Models 

We use two different types of physical models: Bonnor- 
Ebert spheres (Bonnor 1956, Ebert 1955); and Plummer- 
like models (Whitworth & Ward-Thompson 2001). Both 
these models allow an initial density distribution that is 
nearly constant for small radii, but approaches a power law 
at large radii. The Bonnor-Ebert sphere provides a good 
fit to the dust continuum emission of the three PPCs that 
we are modeling (Evans et al. 2001). The Plummer-like 
model can fit the dust continuum data for L1544 (Whit- 
worth & Ward-Thompson 2001). It has the advantage 
that it provides a simple solution for the density and ve- 
locity field during a free-fall collapse from the initial state, 



thereby providing a velocity field, which we need for mod- 
eling some of our lines. The initial Plummer-like radial 
density profile (Whitworth and Ward-Thompson 2001) is 



p(r, t = 0)= Pflat 
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Even though this model is very simple, in that it as- 
sumes a pressure-free collapse, it may be able to explain 
the later phase of the evolution of PPCs. Li et al. (2002) 
show that the evolution of a core driven by ambipolar dif- 
fusion becomes closer to the evolution of a non-magnetic, 
free-falling core at higher density. 

We use the Bonnor-Ebert sphere model that best fit 
the dust emission (Evans et al. 2001), and we assume a 
purely microturbulent velocity field, except for the HCO+ 
J = 3 - 2 lines in L1544 and L1689B and the Hi3co+ 
lines in LI 544. Because those line profiles show deep self- 
absorption features indicative of systematic radial veloc- 
ities, we use a Plummer-like model for them. To con- 
strain the various parameters in the Plummer-like model, 
we modeled the dust emission, as was done by Evans et 
al. (2001). The resulting values that fit best the data for 
both L1689B and L1544 are as follows: 77 = 3, the ini- 
tial pjiaL (t=0) of 5 X 10^ cm-3, and Rfiat of 2800 AU. 
Whitworth and Ward-Thompson found a steeper density 
structure (77 = 4) and a bigger Rfiat (5350 AU) because 
they used an isothermal sphere for their model. When 
the temperature decline toward the center is included, a 
smaller Rfiat is needed. The evolutionary timescale that 
fits the dust emission well is half the timescale at which 
the central point-mass is formed (the start of the Class 
stage). The maximum velocity in this model is about 0.08 
km/s, which is consistent with the results of Tafalla et al. 
(1998) and Williams et al. (1999) in L1544. The resulting 
dcmsity structure is compared to that of the best fitting 
Bonnor-Ebert sphere in Figure 4a; they are very similar. 
Likewise, the resulting dust temperature structure is very 
similar (Figure 4b) . The velocity field from the Plummer- 
like model is shown in Figure 4c. 

While the density and temperature structure of the best- 
fitting Plummer-like model is very similar to that of the 
best fitting Bonnor-Ebert sphere, the former does contain 
a systematic velocity field. To check whether that field af- 
fected the results for lines modeled with the Bonnor-Ebert 
sphere, we modeled all other lines using the Plummer-like 
model. The velocity structure does not affect the line pro- 
files of Ci^O, DCO+, and Hi3C0+; the two models have 
given the same results. 

5.3. Results of Models 

The main purpose of this modeling is to find how much 
a specific molecule is depleted, and how the depletion is 
distributed with radius by simulating molecular line pro- 
files at multiple positions within a core and by comparing 
the simulated line profiles with observed line profiles. We 
tried several functional forms, such as exponential, power 
laws, Plummer-like, and step functions, for the distribu- 
tion of a molecular abundance. In addition, we applied 
the results of Li et al.'s dynamical and chemical model 
(2002) to our modeling, but their model has too high a 
systematic velocity to fit the observed molecular line pro- 
files. However, the group is improving their chemical and 
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dynamical models while communicating with us. Figure 
5 and Table 6 shows the comparison between results of 
different functional forms of depletion with radius in the 
C^^O 2—1 line simulations. We used integrated tempera- 
tures to calculate in each model. Among those distribu- 
tions, a step functional distribution fits best the observed 
line profiles with radius. This seems reasonable because 
the results of chemical models show a sharp decrease in 
abundance with density or time. Here, we mainly show 
the results of the modeling with the abundance distribu- 
tion of a step function that has three free parameters: the 
undepleted abundance of a molecule (Xq), the fractional 
depletion { fo = X/Xq), and the radius inside which a 
molecule is depleted (r^). Each can be different in differ- 
ent species. Therefore, the total depiction of a molecule is 
dependent on the combination of fn and ro- We empha- 
size that the derived depletion factor from this modeling 
has an uncertainty of about a factor of 3 because of the 
uncertainty of the density structure calculated from the 
dust analysis. 

We compared the distributions of observed and mod- 
eled integrated temperatures with radius to find the best 
fits. Figure 6. a shows the distribution of the reduced 
of the models in the C^^O 2-1 and 3-2 lines of L1512 in 
the space of /d and ro- The depletion radius is well con- 
strained around 0.075 pc, but the fractional depletion is 
constrained only to be greater than about 25. This trend 
occurs in every molecule; that is significantly depleted, so 
we consider fu to be the lower limit of the depletion fac- 
tor. In other words, the data are consistent with complete 
depletion inside some rjj. However, the distribution of 
constrains only the upper limit of the depletion radius in 
H13C0+ and DCO+ of L1689B and in DCO+ of L1544, 
which show small r/j and fo (see Table 8 and Table 9). 
For example, the in H^^CO+ of L1689B changes within 
a factor of about 2 (from 3 to 7) in all range of fo within 
rn < 0.011 (Figure 6.b). As we mentioned above, the 
density structure has uncertainty of a factor of about 3, 
so the calculated depletion in this model could have the 
same uncertainty. 

5.3.1. Ci*0 J = 2- lonrfJ = 3- 2 

The best fits of C^^O J = 2- landJ = 3- 2 line 
profiles in the three cores are shown in Figure 7. The pa- 
rameters of the best fit models arc summarized in Tabic 7. 
We used the standard C^'^O abundance of 4.82 x lO^"" for 
the undepleted abundance (Xq) for all cores, which was 
used to calculate the column density of H2 (§ 4.2.1). Ac- 
cording to the results, CO molecules are depleted within 
0.075 pc (about 110") by a factor of 25 in L1512, so CO 
is depleted everywhere in our C^'^O map. In the case of 
L1544, the depletion factor is 25 within 0.045 pc (about 
70"). Ci*0 J = 2- landJ = 3- 2 line profiles simu- 
lated simultaneously at different positions from the cen- 
ter fit the observed line profiles well in L1512 and L1544. 
On the other hand, we could not find a model to fit well 
all the observed line profiles from the center to the off- 
position of 90" in L1689B. The antenna temperatures of 
observed J = 2 — 1 and J = 3 — 2 lines are almost the 
same or increase outward from the center in L1689B, even 
though the outermost lines are not well sampled. This 
might indicate that the large molecular cloud surround- 



ing this source (Loren 1989) contributes to the observed 
C-^^O lines. We tested a core that has a warm envelope 
surrounding an inner, denser core. The core has the same 
density and temperature structures inside 0.15 pc as those 
calculated from dust emission, and a constant density and 
temperature of lO'^ cm~"^ and 50 K from 0.15 to 0.6 pc, 
which are based on the ^^CO observations (Loren 1989). 
The temperature of 50 K is the upper limit in the p Oph 
molecular complex. The modeled line profiles are shown in 
Figure 7c, which fit better the observed line profiles than 
those in the previous model that does not have a warm 
envelope. In this model, the abundance distribution is the 
same as the previous one. In order to avoid confusion by 
the envelope, we need to use a more rare isotope or higher 
J transition to trace the high density core. Even though 
L1512 and L1544 are believed to be well isolated, they 
might also be surrounded by bigger clouds. The modeled 
results, however, fit well the observed C^^O J = 2 — 1 to 
the off-position of 90". The influence of bigger structures 
in L1512 and L1544 does not seem significant compared 
to L1689B. 

Jessop and Ward-Thompson (2001) modeled the C^^O 
J=2— 1 and .1=3—2 lines of L1689B using a density struc- 
ture with the same central density {pfiat) and the same 
Rfiat as calculated by Andre et al. (1996) and a constant 
temperature. After comparing to observations they argued 
that the flat distribution of C^^O is caused by the signifi- 
cant (about 95%) depletion of CO. They pointed out that 
the gas temperature must be above 14 K to have better 
agreement with observations. However, they compared the 
model to the observation within 32", which cannot give a 
good constraint on the outer envelope. In addition, the 
lower density and higher temperature at the inner region 
than those assumed in this study make the transitions op- 
tically thinner than the actual lines, so the depletion factor 
to fit the observations could be overestimated. We mod- 
eled C^^O lines using our density structure and a constant 
temperature of 14 K. The best-fit model has ro ~ 0.035 
and Jd 10, which are bigger than our best-fit model, 
where the temperature structure is calculated from the 
analysis of dust emission. The modeled C-'^^O J=2— 1 lines 
within 60" are stronger than the observed lines. However, 
the line at 90" is still weaker than the observed line even 
though it is stronger than that in our best-fit model with- 
out a warm envelope. 

5.3.2. H'^^CO+ J=l-0andJ = 3-2 

Except for the C^^O lines, Xq is not fixed in model- 
ing molecular lines, but it can be easily decided by the 
outermost line profile. The results of the best models in 
H13C0+ J = 1 - and J = 3 - 2 are shown in Figure 8 
and Table 8. We have observed the Hi3C0+ J = 1 - line 
with a 20" grid out to the distance of 60" from the center 
with an 18" beam. The J = 3 — 2 line has been observed 
with 30" spacing with a 26" beam. In L1512, the Hi^''CO+ 
J = 3 — 2 line, whose critical density is about 10^ cm~^, 
was not detected, so we simulated only the J = 1 — line. 
The depletion factor in L1512 is 25 within 0.021 pc (38"). 
In L1689B, the best-fit model of J = 1 — line also fits 
J = 3 — 2 line very well. The depletion factor and radius 
in L1689B are much smaller than those in L1512. How- 
ever, the undepleted abundance of H^^CO+ in L1512 is 
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1.4 times bigger than in L1689B. 

The observed H^^CO+ J = 1 — hnc, wliose critical 
density is about 10^ cm~^, shows a blue skewed profile in 
LI 544. Contracting cores show stronger blue peaks and 
self-absorption dips in optically thick lines (Myers ct al. 
1996). Evidence of inward motion in L1544 has been ob- 
served by Tafalla et al. (1998) and Williams et al. (1999). 
Therefore, we used the best-fit Plummer-like model in the 
simulation of the H^^CO+ J = 1-0 line. In L1544, models 
fit well the shapes of the observed H^'^CO+ J = 1 — line 
profiles, but the modeled J = 3 — 2 line is s^ron^er than the 
observed line by a factor of 2. The rms noise of this line 
is about 0.015 K {5a detection), and the baseline is poor, 
so the discrepancy may not be very significant. There are 
several possibilities. If there is molecular gas surrounding 
LI 544, this less dense gas would increase the J = 1 — line 
more than the J = 3 — 2 line. Therefore, our models would 
overestimate the abundance of the inner, denser core to fit 
the J = 1 — line. As a result, the modeled J = 3 — 2 line 
can be stronger than the observed line. Another possible 
explanation is that the depletion fraction increases toward 
the denser, inner region. Since the J = 3 — 2 line traces 
denser gas than the J = 1 — line, the J = 3 — 2 line 
in our models where we used a step function of the abun- 
dance would be stronger than the observed line if H^^CO+ 
is depleted more at the center. The uncertainties in the 
observed line or in the density and temperature structures 
could explain the difference between the modeled and ob- 
served lines. 

5.3.3. DCO+ J = 3-2 

The modeling of DCO+ J = 3 — 2 does not constrain 
free parameters well because the observations in the line 
do not cover positions far from the center. The results 
are shown in Figure 9 and in Table 9. In L1512, only the 
center position was observed, so we just used a constant 
abundance to fit the line profile. The abundance of DCO+ 
that fits the line profile well is 2.8 x 10"^'^. However, the 
line profiles of two positions in LI 544 and three positions 
in L1689B were compared with the results of models. The 
factor and the radius of the depletion of DCO+ molecule 
are much smaller than those of CO molecules calculated 
by simulating C^^O J = 2 - 1 and J = 3 - 2 lines. This 
result is consistent with the results of the chemical model 
of Caselli et al. (2002). The radius of the depletion is less 
(20") than the beam size (35") in L1689B and similar to 
the beam size in L1544. 

5.3.4. HCO+ J = 3-2 

Infall asymmetry in PPCs has been observed (Tafalla et 
al. 1998, Lee et al. 1999, Gregersen and Evans 2000) in op- 
tically thick lines. If the gas of the outermost envelope of 
a core that has inward motion has sufficient optical depth 
and has a lower excitation temperature than the gas of 
the inner envelope in a given molecular line, the molecular 
line shows infall asymmetry, in which the intensity ratio 
of blue peak to red peak increases with time (Gregersen 
et al. 1997). We model one of the optically thick lines, 
HCO+ J = 3 — 2 in the three PPCs and compare the re- 
sults with the observational results of Gregersen and Evans 
(2000). In order to model an HCO+ J = 3 - 2 line that 
has infall asymmetry, we used the best-fit Plummer-like 



model. Figure 10 shows the difference in MC modeling of 
the HCO+ J = 3 — 2 line for two different physical models: 
Bonnor-Ebert sphere and Plummer-like model. The latter 
includes an infall velocity structure and produces a blue 
asymmetry. The HCO+ J = 3 — 2 line in L1512 does not 
show a self-absorption dip and asymmetric profile so that 
we model the line in L1512 using the best-fit Bonnor-Ebert 
sphere. 

The results of modeling HCO"'" J = 3 — 2 in the three 
PPCs are shown in Figure 11, and in Table 10. Because 
L1512 has only one spectrum at the position 30" away from 
the center, we used a constant abundance of HCO+ for this 
core. The abundance of the best fit is about 1 x 10~^. In 
L1544, the intensity of the HCO+ J = 3 - 2 line at 30" is 
half that at the center, unlike the H^'^CO"'" 1-0 line, whose 
intensity is almost flat to 40". As a result, an exponential 
abundance distribution as well as a step function also fits 
well the hne profiles. In L1689B, the HCO+ J = 3 - 2 
line has a stronger red-peak asymmetry at the center and 
15" away, but it has a stronger blue-peak asymmetry at 
30" and 45". This complication of line profiles causes dif- 
ficulty in fitting the exact line profiles in this core. There- 
fore, we tried to fit the integrated temperatures of the lines 
with radius. The exponential function of the abundance 
shows similar results to the step function. In order to test 
the functional form of the abundance, we need a more 
extended map. The interesting result is that the HCO"'" 
J = 3 — 2 line of L1689B has a similar strength to that of 
L1544, unlike H^^CO"^ 1-0, which is stronger in L1689B 
than in L1544. This might be due to the fact that L1689B 
is embedded in a larger molecular cloud that causes more 
absorption of the HCO+ 3—2 line. 

On the whole, the abundance of HCO"*" in this result is 
very small compared to the abundance that is expected 
from H^^CO"*" and the typical ratio of C and ^^C in the 
local interstellar medium of about 77 (Wilson and Rood 
1994). However, our results show the abundance ratio of 
HCO"'" and H^^CO"'" is approximately from 3 to 15 in the 
three PPCs. (We do not think that the ratio in L1512 is 
correct because we used only one spectrum in each line 
to calculate the abimdance.) In addition, the depletion 
factors and the depletion radii calculated from the HCO"*" 
line are different from those calculated from the H^^CO"*" 
lines. These differences and the low abundance ratio may 
result from the gas components outside the cores that we 
modeled in this study because the HCO+ line can trace 
less dense gas than H^^CO"*". As a result, the HCO"'" line 
can be absorbed by the surrounding material significantly, 
but the H^'^CO"'" line may be not absorbed as much as the 
HCO+ line. In our models, we limited the size of every core 
to 0.15 pc so that our models cannot account for the gas 
outside the core. Then, we would have underestimated the 
abundance trying to fit a line that has been weakened by 
the absorption by the surrounding gas component. How- 
ever, this explanation is just one possibility. We cannot 
ignore the uncertainties of many other parameters, so we 
leave this as future work. 

6. DISCUSSION 

6.1. Depletion 

We concluded that CO is significantly depleted in L1512 
and L1544. It is difficult to estimate the depletion of CO 
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in L1689B because of the gas surrounding this core and 
the optical depths of the lines used, but the depletion ap- 
pears to be much less. The HCO+ molecule is depleted in 
L1512 and L1544, but we do not see significant depletion of 
this molecule in L1689B, based on modeling the H^''CO+ 
J = 1 — line. The optical depth and possible surround- 
ing material also affect the calculation of the depletion of 
HCO+ in L1544. DCO+ is not significantly depleted in 
any of the PPCs. This result is consistent with results in 
many other starless cores where the distribution of DCO+ 
fits well the distribution of dust continuum (Myers 2002, 
personal communication). However, we need to use better 
resolution to do a more precise study because the calcu- 
lated depletion radius is not bigger than the actual beam 
size (35"). Caselli et al. (2002) show that the depletion 
radius of DCO+ is about 3000 AU in L1544. 

The analysis of the integrated intensity shows that the 
earliest-time molecule CCS is not peaked toward the cen- 
ters of any of the PPCs, imlike the dust emission. Thus, 
CCS is depleted significantly in all three. As we saw, how- 
ever, the CCS intensity in L1544 is twice stronger than 
that in L1689B even though the dust emission of L1689B 
and L1544 is similar. This result suggests that the deple- 
tion of CCS is connected to the depletion of CO. Li et al. 
(2002) modeled chemical abundance changes in PPCs us- 
ing a dynamical model with a magnetic field, which drives 
the ambipolar diffusion, and they got a relatively high 
abundance of CCS. They indicated that this abundance 
could result from the interplay between significant deple- 
tion of CO molecules and late-time hydrocarbon chem- 
istry. Ruffle et al. (1997, 1999) also showed the same 
result in their models of the time-dependence of several 
molecules. According to their models, the abundance of 
CCH and CCS have later secondary maxima only when 
the freeze-out timescale is long compared to the chemical 
timescale or the collapse timescale. CCS increases quickly 
to the first maximum of its abundance in time and is de- 
pleted very early as the density increases because CCS is 
destroyed by increasing HCO+ or H3O+, and this polar 
molecule can be tightly bound in grain mantles. How- 
ever, in the region where CO molecules are depleted sig- 
nificantly, the C"*" ion reacts more with molecules that are 
not oxygen-bearing species so that the precursors of CCS 
such as CCH increase again. This results in a second in- 
crease of the CCS abundance. The trend has not shown up 
in the models of Bergin & Langer (1997) and Aikawa et al. 
(2001) because the models used dense initial conditions in 
which a core evolves very quickly dynamically, and CCS 
did not show the second maxima. We think that the CCS 
molecule is possibly in the second enhancement in the late 
stage of chemical evolution resulting from the significant 
depletion of CO in L1512 and L1544. The integrated in- 
tensity of the CCS line shows bumps in the center regions 
of L1512 and L1544 (Figure 3). If CCS is depleted with 
time without the second enhancement, we should see holes 
rather than bumps. The other possible explanation of the 
enhancement of CCS is mixing of some atomic carbon from 
a translucent envelope into the core (van Dishoeck 2002, 
personal communication) . 

According to the comparison between the intensity of 
N2H"'" lines and the dust emission intensity, this molecule is 
possibly depleted in L1512, but not in L1544 and L1689B, 



where its intensity follows well the dust emission inten- 
sity. N2H+ is known as a very late-time molecule, but 
Bergin et al (2002) have shown that N2H+ is depleted in a 
cold dark cloud, B68. The weak intensity of N2H+ line in 
L1689B, compared to in L1544, might indicate that this 
molecule has not yet reached its maximum because of less 
abundant precursors (N2 and H3") but a more abundant 
destroyer (CO). 

Bergin and Langer (1997) modeled chemical evolution in 
PPCs including gas-grain interactions as well as gas-phase 
reactions and compared two kinds of grain mantle proper- 
ties: a weakly bound CO mantle and tightly bound H2O 
mantle. Their results show that sulfur-bearing molecules 
such as CS and SO are very sensitive to density and de- 
pleted seriously as the density of a core increases, regard- 
less of the properties of dust grains. On the other hand, 
CO and HCO+ molecules are depleted in the model with 
H2O grain mantles but remain in gas phase in the model 
with CO grain mantles. In contrast, NH3 and N2H"'" do 
not show depletion in any model because of the low bind- 
ing energy of N2 to H2O and CO mantles. Therefore, the 
comparison of the results of their chemical models to our 
results indicates that the dust grains in the three PPCs 
have tightly bound H2O mantles. In addition, Aikawa 
et al. (2001) predicted the abundance distributions of 
molecules by numerical chemical models simply by apply- 
ing a constant delay factor to a Larson-Penston model and 
considering different sets of adsorption energies, and they 
compared their results with the observed molecular abun- 
dances in L1544. According to their results, CCS and CO 
molecules are depleted, but N2H+ and NH3 molecules are 
more abundant in more slowly collapsing cores at the den- 
sity peak. 

To summarize our results on depletion, CCS, CO, and 
HCO"*" are significantly depleted in L1512 and in L1544. 
On the other hand, the depletion of DCO+ in the two cores 
is not significant. N2H+ is not depleted in L1544, but it 
is possibly depleted in L1512. In contrast, only CCS is 
significantly depleted in L1689B. It is hard for us to esti- 
mate the depletion of CO in L1689B because of the optical 
depth eff'ect and the material surrounding the cores. 

6.2. Ionization 

The coupling between ions and the magnetic field is cru- 
cial in regulating star formation rates if cores are magnet- 
ically subcritical. The timescale of ambipolar diffusion de- 
pends on the ionization fraction (tAD ~ 2.5x 10"'^^a;(e), Shu 
et al. 1987). In addition, the gas-phase chemical reactions 
are also dependent on the ionization fraction. Therefore, 
calculating the ionization fraction accurately is significant 
for understanding the chemistry as well as the timescale 
for ambipolar diffusion. However, it is not easy to calcu- 
late the ionization fraction in protostellar cores because it 
depends on the depletions of molecules as well as density 
structures and cosmic-ray ionization rate. 

We used equation (16) in Caselli (2002) to calculate the 
ionization fractions: 
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where fo is the depletion factor of CO, x(0) is the abun- 
dance of atomic oxygen (1.5 x 10~^, Caselli et al. 1998), 
f^ is the depletion factor of O, and [DCO+]/[HCO+] is the 
abundance ratio of two molecules, which can be calculated 
by comparing the column densities of DCO+ and H^^CO"'" 
and using 77 (±7) for the ^^C/^^C ratio (Wilson & Rood 
1994). Caselli et al. (1998) and Caselh (2002) derived this 
equation assuming that (i) the main source of HCO"*" is 
the reaction between and CO, (ii) the deuterium frac- 
tionation is due to the reaction between H3 and HD, and 
(iii) molecular ions are destroyed mainly by electrons and 
neutral species such as CO and O. In fact, Caselli et al. 
(2002) have showed that N2H+ and N2D+ are better trac- 
ers of ionization fraction because HCO"'" and DCO"*" are 
possibly depicted in the central region. However, we do 
not have N2D+ data, so we use DCO"'" and HCO+, which 
can give correct information in the region with r > 3000 
AU. 

We calculate, for each line of sight, the average x{e) 
along the line of sight. This approach is not completely 
self-consistent because the method is only valid for homo- 
geneous clouds. We use our model to simulate observa- 
tion with the same beam size. The column density of each 
molecule has been calculated from the simulated line based 
on the assumption of an optical thin line. Figure 12 shows 
the ionization fraction with radius in the three PPCs. 

There are some caveats in this calculation. First, we 
could not assess the depletion factor of CO in L1689B cor- 
rectly because of the optical depth problem and the pos- 
sible material surrounding the core (§ 5.2.1). Second, the 
model of H^^CO"'" lines in L1544 does not match the ob- 
servations very well so that the calculated column density 
from the simulated H^^CO"'" J = 3 — 2 has more uncer- 
tainty. Third, since we do not have enough observed data 
to assess the depletion radius of DCO"*" in L1512, the col- 
umn density of DCO"*" is uncertain in this core too. In 
addition to the caveats related to models, the equation 
that we used to calculate ionization fraction includes the 
abundance and the depiction fraction of O that arc not 
well known. We supposed that the depletion of atomic 
oxygen might be proportional to the depletion factors of 
other molecules. Wc used f^ of 3, 2, and 1 for L1512, 
L1544, and L1689B, respectively. 

In spite of those caveats, the three cores have similar 
average ionization fractions (about 5 x lO""") along the 
line of sight toward the centers. These results are similar 
to the ionization fractions calculated in other protostellar 
cores (Caselli et al. 1998) but much bigger than the result 
(about 10"" at the center) of Caselli et al. (2002). Caselli 
et al. (1998) have used the integrated temperatures of 
the actual line profiles, but Caselli et al. (2002) have used 
the results of their "L1544 best-fit" chemical model, which 
takes into account the cloud density structure, molecu- 
lar freeze out, and the recombination of molecular ions 
on grain surfaces, in order to infer the electron fraction 
as a function of cloud radius. Our data do not rule out 
still smaller values of x{e) at small radii. The ionization 
fraction in L1512 and L1544 decreases toward the center, 
which trend has been shown in Caselli et al. (2002) and 
Caselli (2002). In contrast, in L1689B, it increases toward 
the center. As a result, the ionization fraction of L1689B 
around 0.06 pc is smaller than those of L1512 and L1544 



by a factor of 10. 

6.3. Timescales 

The timcscalc of ambipolar diffusion (tAD ~ 2.5 x 
10^'^a;(e), Shu et al. 1987) calculated from the above ion- 
ization fraction toward the centers is about 1.25 x 10^ 
years. This timescale applies to cores that are magneti- 
cally subcritical. This ambipolar timescale represents the 
duration of contraction from the present stage to the stage 
of Class 0. These are much bigger than the free fall 
timescale (« 10^ years). However, x{e) could be much 
smaller at the center. Even though these three cores have 
similar ambipolar diffusion timescale toward the centers, 
our results suggest that they have different chemical status 
from each other. 

L1512 seems to have evolved chemically more than 
L1544 showing that CCS, CO, HCO+, and possiblely 
N2H+ are depleted. Caselli (2002) modeled a core that 
is less dense and less centrally concentrated than L1544 
to show the chemistry in an early stage in the dynami- 
cal evolution. According to their result, molecules in the 
model core are not significantly depleted. The model core 
is similar to L1512 in density structure, but surprisingly, 
L1512 shows a bigger depletion radius of CO than L1544. 
In addition, N2H+ at the center of L1512 might start to 
be depleted. This result could result from the slower con- 
traction of L1512 than L1544 so that L1512 could have a 
longer time to evolve chemically. Another possibility is a 
more eflacient molecular freeze out in L1512 because the 
dust grains in this object are smaller than in the other two 
cores, so that the depletion timescale in L1512 is shorter 
(van Dishoeck et al. 1993; Caselli et al., in prep.). This 
could be possible because L1512 is less dense than the 
other two cores and has not have enough time for grain 
coagulation (Evans et al. 2001). 

On the other hand, L1689B has not evolved much chem- 
ically even though the density structure is very similar to 
LI 544 so that it is in a later stage of dynamical evolution 
than is L1512. This tells us that the chemical evolution 
is not a simple function of density, unlike the results of 
other chemical models. According to this result, we should 
not ignore the possibility that L1689B could be magnet- 
ically supercritical, and the core has evolved faster than 
the other two sources dynamically so that L1689B has not 
had enough time to evolve chemically. The broader lines in 
L1689B (Figure 1 and § 3.1) indicate more active kinemat- 
ics than the other two sources. The total masses of L1544 
and L1689B are similar according to the dust modeling. 
Therefore, L1689B has to have a weaker magnetic field or 
a smaller ionization fraction in order to be magnetically 
supercritical. Crutcher & Troland (2000) measured the 
magnetic field along the fine of sight in L1544 (i?LOS ~ H 
/xG) by using the Zeeman effect of the OH lines. There 
is no specific measurement of magnetic fields in L1689B. 
However, Troland et al. (1996) measured a OH Zeeman ef- 
fect in the p Oph cloud indicating i? < 10 fiG. In addition, 
our calculation shows the ionization fraction in L1689B is 
smaller than in the other two cores in the outer part, where 
the calculation is more credible. The other possibility of 
the faster dynamical evolution of L1689B than L1544 is a 
larger external pressure (Galli et al. 2002) because L1689B 
is in the p Oph complex (see §5.3.1). 
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Based on our results, we suggest that the stage of dy- 
namical evolution of a core cannot be simply probed by its 
chemical status because the chemical evolution depends 
on the size of dust grains or the absolute timescale, during 
which the core has been in a given environment, as well 
as its density structure that shows the relative stage of 
dynamical evolution. We summarize the dynamical and 
chemical evolutionary stages of each core in Table 11. 

7. SUMMARY 

The summary of our results is as follows: 

1) The difference between the column densities inferred 
from dust emission and the omission of a given molecu- 
lar line can partially result from the optical depth of the 
line. However, the depletion of a molecule must explain 
the remaining difference. In L1689B, the main cause of the 
difference between the dust emission and the CO emission 
is the optical depth in the lines. However, in L1512, the 
depletion of CO is the major source of the difference be- 
tween the dust emission and the CO emission. 

2) In our Monte Carlo simulation of molecular lines, a step 
function with depletion in the inner region produces the 
best fit to the abundance variation of a molecule in PPCs. 

3) The CO and HCO"*" molecules are depleted significantly 
in L1512 and L1544, but not in L1689B. The gas around 
L1689B makes it difficult to calculate the actual depletion 
of CO. 

4) The depletion of DCO+ is not substantial in the three 
PPCs. 

5) The CCS molecule is substantially depleted in the three 
PPCs. However, CCS might be in the second enhancement 
in L1512 and LI 544 because of the significant depletion of 
CO. 

6) The distribution of the integrated intensity of N2H+ 
J = 1 — suggests the possible depletion of N2H+ in L1512. 
On the other hand, N2H+ follows dust emission very well 
in L1544. In L1689B, N2H+ has not had enough time to 
reach its maximum. 

7) The velocity structures calculated by a simple dynam- 



ical model that uses a Phimmer-like density profile and 
a free-fall collapse model can produce the asymmetry of 
HCO+ J = 3 - 2 fines in L1544 and L1689B. However, 
the abundance calculated by modeling the lines is much 
smaller than what is expected from modeling H-'^'^CO^. 
The surrounding gas could possibly account for the low 
HCO+ abundance of the models. 

8) The ionization fraction is similar toward the centers of 
the three cores even though L1512 is less centrally con- 
densed than L1544 and L1689B. However, L1512 has a 
bigger depletion radius of CO than L1544 and possible 
depletion of N2H+ at the center. In contrast, L1689B 
does not show the maturity in chemical evolution com- 
pared to the other two cores. This difference suggests that 
the chemical evolution depends on the size of dust grains 
or the absolute timescale during which a core stays in a 
given environment as well as the density structure of the 
core that shows the relative dynamical evolutionary stage. 

9) The evolved density structure and the young chemistry 
of L1689B suggest that this core is possibly experiencing 
a free-fall collapse rather than ambipolar diffusion. As 
a result, L1689B evolves dynamically too fast to evolve 
chemically. This would be possible if the magnetic field 
in L1689B is weak enough, or the external pressure is big 
enough to make L1689B magnetically supercritical. 
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Fig. 1. — The observed spectra toward the dust peaks of three PPCs. The dotted line in the boxs of the third row shows the results of 
fitting the hyperfine structure of C^'^O J = 2 — 1 line profiles using solutions in CLASS. 



14 



100 



50 



■50 



100 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 

dotted: C180 2-1 L1512 



0.50 



0.26 





0.58 



0.56 



00 



50 





Aa 



■50 



100 



Fic;. 2. A. — These maps show the distributions of 850 fim dust continuum intensity (grey scale, Shirley et al. 2000) and integrated intensities 
of C-'^^O J = 2 — 1 line (dotted contour scale) and DCO+ J = 3 — 2 line (solid contour scale). The (0,0) position is the reference position in 
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J = 2 - 1 and &''0 J = 2 - 1 transitions. 



15 




16 




17 



I I I I I I I I I I I I I I I I I I I I I I 
L1512 



I ' ' ' ' I 



O 

J -1 




I I I I I I I I I I I I I I I I I I 



LI 544 



I I I I I I 




H13C0+ 



I I I I I I I I I 



L1689B 




I I I I I I I I I I I I I I I I I I 



22.5 



22 



o 



21.5 



21 






20.5 - 



■ N(H,),,„ 

N(Hg)QjgQ 



_L 



_L 



_L 



_L 



100 50 



-50 -100 100 50 -50 -100 100 50 
a offset (arcsec) 



-50 -100 



Fig. 3. — These plots in upper panels compare the integrated intensities of CCS Nj = 43 — 82, H^^CO+ J = 1 — 0, and N2H+ J = 1 — 
lines with Ssbo through the cuts marked in Figure 2. Here, Ss50 is shifted by 0.9, and /hi^co-i- shifted by 0.4 in the logarithmic scale. 
These plots in lower panels compare the column densities of hydrogen molecules calculated from the 850 /im dust continuum emission, C^®0 
J = 2 — 1 line (without the correction of t), and C'^^O J = 2 — 1 line (with correction of r) through the cuts marked in Figure 2. The errors 
in molecular intensities only account for the rms noises of spectra so that the error bars show the minimum errors. In N(H2) derived from 
C^'^O of L1544 and L1689B, the error in fitting the hyperfine structure is also included. 



18 




Fig. 4. — The comparison between physical parameters of the Bonnor-Ebert model and the Plummer-like model; (a) the density distribution, 
(b) the temperature distribution, and (c) the velocity distribution. In the density structure of Plummer-like model, we forced the lower limit 
of density to be 10'^ cm~^. 
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Fig. 6.A. — The distribution of the reduced of models that have different ro and fo in C^*0 2—1 and 3—2 toward L1512. The contour 
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Fig. 6.B.— The distribution of the reduced of models that have different td and fo in Hi3cO+ 1-0 and 3-2 toward L1689B. The 
contour levels are from 3 to 33 increasing by the interval of 2. 
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Fig. 7.B.— The results of MC modeling in C^^O J = 3 - 2. 
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Fig. 7.C.— The results of MC modeling of L1689B in C^*0 lines including a warm envelope. 
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Fig. 8.A. — The results of MC modeling in 11^300+ J = 1 - 0. Hi3CO+ J = 1 - and J = 3 - 2 lines have been modeled simultaneously, 
and we chose the abundance structure that fitted both of the lines. In order to model these line profiles, in L1512 and L1689B, we used the 
density and temperature profiles given by Bonnor-Ebert dust models (Evans et al. 2001), but we calculated the structures of density and 
temperature in L1544 combining the empirical dynamics model of Whitworth and Ward-Thompson (2001) and dust radiative transfer code. 
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Fig. 8.B.— The results of MC modeling in Hi3cO+ J = 3 - 2. 
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Fig. 9. — The results of MC modeling in DCO+ J = 3 - 2 line. 
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Fig. 10. — The comparison of MC modeling in HCO+ .1 = 3 — 2 line of L1544 for two different physical models: Bonnor-Ebert sphere 
and Plummcr-likc model (sec Figure 4). Both have the same abundance structure, but only the Plummer- like model, with an infall velocity 
structure, produces the blue asymmetry shown in observed line profiles. 
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Fig. 11. — The results of MC modeling in HCO+ J = 3 — 2 line. In L1544 and L1689B, the upper and the lower panels show the results 
by the step function and the exponential function of HCO+ abundance, respectively. In L1512, we modeled the line profile using a constant 
abundance. 
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Fig. 12. — The ionization fractions in three PPCs. The error bars represent the amount of uncertainty in the ionization fraction when the 
ratio of DCO+/HCO+ and the depletion fraction of CO molecule changes up to 30 % of the calculated values in each point. 
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Table 1 
List of sources 



Name 


R.A. (1950.0) 
hh mm ss 


Dec. (1950.0) 

o / // 


Vlsr 
km 


distance 
pc 


offset of dust peak 
arcsec 


L1512 


05 00 54.4 


32 39 37 


7.09 


140 


(-19, -26) 


L1544 


05 01 13.1 


25 06 36 


7.12 


140 


(-1, 4) 


L1689B 


16 31 47.0 


-24 57 22 


3.61 


125 


(-1, -12) 



Table 2 
List of observed lines 



Molecule 


Transition 


Frequency 


Beam Width 


Velocity Resolution 








MHz 


arcsec 


km s~^ 


L1512/L1544/L1689B 




J = 2-1"* 


224714.385 


33 


0.17 


0.65/0.65/0.80 




J = 2-P 


219560.357 


34 


0.15 


0.69/0.66/0.60 




J = 3 - 2^ 


329330.507 


26 


0.10 


0.59/0.59/0.63 


HCO+ 


J = 3 - 2^'= 


267557.620 


26 


0.15 


0.49/0.66/0.65 


H13C0+ 


J = 1 - 0'' 


86754.330 


18 


0.14 


0.49/0.49/0.49 




J = 3 - 2^'^ 


260255.617 


26 


0.16 


none/0.49/0.65 


DCO+ 


J = 3 - 2^ 


216112.605 


35 


0.15 


0.65/0.69/0.60 


CCS 


Nj = 43 - 32'' 


45379.033 


37 


0.27 


0.77/0.77/0.77 


N2H+ 


J = 1 - O'' 


93173.809 


17 


0.13 


0.50/0.50/0.50 



"^observed with CSO 10.4 m telescope 
''observed with Nobeyama 45 m telescope 
'^Gregersen and Evans, 2000 

'^'qmb = T\/T^, where is the radiation temperature, and Tr = only if a source fills the beam 
(Kutner and Ulich 1981). 



Table 3 
Parameters for molecules 



Isotope 


Rotational constant 


Dipole moment 


Abundance 




B (MHz) 


(Debye) 


X{= rii/nH^) 


CiSQ 


54891.423 


0.11079 


4.8 X 10"^ 


C"0 


56179.982 


0.11034 


1.5 X 10"^ 



Table 4 

Optical depth of C^^'O J = 2 — 1 line at center position 



Source 


Total optical depth 


L1512 


0.10 (±2.17) 


L1544 


1.03 (±0.42) 


L1689B 


0.79 (±0.24) 
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Table 5 

Depletion factor of CO molecule at center position 



Source N{ll2)Ss,o/Ni^2)ciro^ 

L1512 11 

L1544 9 
L1689B 2 



'^We have corrected for the optical 
depth effect to calculate the column 
density. 



Table 6 

The in modeling the C^*0 2—1 lines of L1544 by the different functional forms of depletion 



Function 



X 



Step 



Plummer-like Xi 

Power-law 
Exponential 



i/r-2 + (0.05)2 

Xoexpi-^) 



X 



Xo'' (r > 0.045 pc) 10 
Xo/25 (r < 0.045 pc) 



76 

117 
93 



"■Xo = 4.82 X 10"^ 



Table 7 

The parameters of best-fit models in C^^O lines 



Source 










^^icroturbulence 








pc 




km S""'^ 


L1512 


4.82 X 10- 


-7 


0.075 


25 


0.15 


L1544 


4.82 X 10- 


-1 


0.045 


25 


0.17 


L1689B'' 


4.82 X 10" 


-7 


0.03 


3 


0.20 



'^X = Xo (r > Tu) and X = Xq/ fo (r < td) where 
Xq is undepleted abundance, is the radius within 
which a molecule is depleted, and fo is the fraction of 
depletion. 

^We could not find a good model to fit all line profiles 
in all observed positions in L1689B. 



Table 8 

The parameters of best-fit models in H^^C0+ lines 



Source 


^0 






fD 


^^icroturbulence 








pc 




km s"^ 


L1512 


3.0 X 10- 


10 


0.021 


25 


0.07 


L1544 


5.0 X 10- 


10 


0.026 


20 


0.20 


L1689B 


2.2 X 10- 


10 


0.012 


5 


0.20 
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Table 9 

The parameters of best-fit models in DCO"*" lines 



Source 


^0 






fD 


^^icroturbulence 








pc 




km S"^ 


L1512 


2.8 X 10" 


-10 






0.09 


L1544 


3.5 X 10- 


-10 


0.022 


3 


0.12 


L1689B 


3.5 X 10- 


-10 


0.011 


4 


0.20 



Table 10 

The I'AHAAIETEHS OF LiEST-FlT MODELS L\ HCO LINES 



Source 




step function 




exponential function*^ 




^^nicroturbulence 




^0 




fD 


^0 


td 








pc 






pc 


km 


L1512 


1.0 X 10"^ 










0.1 


L1544 


3.0 X 10-*^ 


0.028 


10 


6.0 X 10^3 


0.04 


0.15 


L1689B 


3.2 X 10-3 


0.023 


15 


6.5 X 10-9 


0.03 


0.20 



X = Xo expi-^) 



Table 11 
The results of this study 





L1512 


L1544 


L1689B 


Degree of central condensation 


Young 


Evolved 


Evolved 


Chemical evolution (depletion of molecules) 


Evolved 


Evolved 


Young 


Timcscalc 


'^che '^dyn 


'Tche ~ '^dy7i 


Tche ^ ^ "^dyn 


Dynamical state 


Stable (?) 


AD (?) 


Free-fall (?) 



